0.1 Import Data

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, broom, skimr, devtools, ggpubr, mgcv, extrafont, mgcViz, here)

theme_set(theme_pubclean(base_size = 14))

# load the function to print GAM figures
source(here("R", "p_gam.R"))
# import
lesion <- read_csv(here("cache", "summarised_lesion_data.csv"))
weather <- read_csv(here("cache", "weather_summary.csv"))

dat <- left_join(lesion, weather, by = c("site", "rep"))

0.2 Fit GAMs

For reproducibility purposes, use set.seed().

set.seed(27)

0.2.1 mod1 - s(Distance)

mod1 <-
   gam(
      mean_pot_count ~ s(distance, k = 5),
      data = dat
   )

summary(mod1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   1.0802     0.0475    22.7 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F             p-value    
## s(distance) 3.93      4 78.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.482   Deviance explained = 48.8%
## GCV = 0.76522  Scale est. = 0.75394   n = 334
print(p_gam(x = getViz(mod1)) +
         ggtitle("s(Distance)"),
      pages = 1)

0.2.2 mod2 - s(Distance) + Precipitation

mod2 <-
   gam(
      mean_pot_count ~ sum_rain+ s(distance, k = 5),
      data = dat
   )

summary(mod2)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ sum_rain + s(distance, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  1.12495    0.05569   20.20 <0.0000000000000002 ***
## sum_rain    -0.01085    0.00709   -1.53                0.13    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F             p-value    
## s(distance) 3.93      4 78.7 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.484   Deviance explained = 49.1%
## GCV = 0.76443  Scale est. = 0.75086   n = 334
print(p_gam(x = getViz(mod2)) +
         ggtitle("s(Distance) + Precipitation"),
      pages = 1)

0.2.3 mod3 - s(Distance) + Windspeed

mod3 <-
   gam(mean_pot_count ~ mws + s(distance, k = 5),
       data = dat)

summary(mod3)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ mws + s(distance, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   0.6440     0.1182    5.45 0.0000001 ***
## mws           0.1227     0.0306    4.01 0.0000747 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df  F             p-value    
## s(distance) 3.93      4 82 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.504   Deviance explained = 51.2%
## GCV = 0.73389  Scale est. = 0.72086   n = 334
print(p_gam(x = getViz(mod3)) +
         ggtitle("s(Distance) + Windspeed"),
      pages = 1)

0.2.4 mod4 - s(Distance) + Windspeed + Precipitation

mod4 <-
   gam(mean_pot_count ~ sum_rain + mws + s(distance, k = 5),
       data = dat)

summary(mod4)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ sum_rain + mws + s(distance, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)  0.67473    0.11852    5.69 0.000000028 ***
## sum_rain    -0.01472    0.00697   -2.11       0.035 *  
## mws          0.13115    0.03069    4.27 0.000025325 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F             p-value    
## s(distance) 3.93      4 82.8 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.509   Deviance explained = 51.8%
## GCV = 0.72844  Scale est. = 0.71333   n = 334
print(p_gam(x = getViz(mod4)) +
         ggtitle("s(Distance) + Windspeed + Precipitation"),
      pages = 1)

0.2.5 mod5 - s(Distance + Windspeed) + Precipitation

mod5 <-
   gam(
      mean_pot_count ~ sum_rain + s(distance + mws, k = 5),
      data = dat
   )
## Warning in term[i] <- attr(terms(reformulate(term[i])), "term.labels"): number
## of items to replace is not a multiple of replacement length
summary(mod5)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ sum_rain + s(distance + mws, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)  1.12495    0.05569   20.20 <0.0000000000000002 ***
## sum_rain    -0.01085    0.00709   -1.53                0.13    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F             p-value    
## s(distance) 3.93      4 78.7 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.484   Deviance explained = 49.1%
## GCV = 0.76443  Scale est. = 0.75086   n = 334
print(p_gam(x = getViz(mod5)) +
         ggtitle("s(Distance + Windspeed) + Precipitation"),
      pages = 1)

0.2.6 mod6 - s(Distance) + s(Windspeed) + Precipitation

mod6 <-
   gam(
      mean_pot_count ~ sum_rain + s(distance, k = 5) + s(mws, k = 5),
      data = dat
   )

summary(mod6)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ sum_rain + s(distance, k = 5) + s(mws, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   1.1879     0.0715   16.61 <0.0000000000000002 ***
## sum_rain     -0.0261     0.0138   -1.89               0.059 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F              p-value    
## s(distance) 3.94      4 93.8 < 0.0000000000000002 ***
## s(mws)      3.93      4 16.0       0.000000000003 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.566   Deviance explained = 57.8%
## GCV = 0.64971  Scale est. = 0.63051   n = 334
print(p_gam(x = getViz(mod6)) +
         ggtitle("s(Distance) + s(Windspeed) + Precipitation"),
      pages = 1)

0.2.7 mod7 - s(Distance) + s(Windspeed)

mod7 <-
   gam(
      mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5),
      data = dat
   )

summary(mod7)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   1.0802     0.0437    24.7 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df  F              p-value    
## s(distance) 3.94   4.00 93 < 0.0000000000000002 ***
## s(mws)      3.92   3.99 16       0.000000000006 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.562   Deviance explained = 57.3%
## GCV = 0.65392  Scale est. = 0.63659   n = 334
print(p_gam(x = getViz(mod7)) +
         ggtitle("s(Distance) + s(Windspeed)"),
      pages = 1)

0.2.8 mod8 - s(Distance) + s(Windspeed) + s(Precipitation)

mod8 <-
   gam(
      mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5) + s(sum_rain, k = 5),
      data = dat
   )

summary(mod8)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5) + s(sum_rain, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   1.0802     0.0434    24.9 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F             p-value    
## s(distance) 3.94   4.00 93.80 <0.0000000000000002 ***
## s(mws)      2.98   3.05  6.71               0.016 *  
## s(sum_rain) 1.93   1.94  1.98               0.190    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.566   Deviance explained = 57.8%
## GCV = 0.64966  Scale est. = 0.63051   n = 334
print(p_gam(x = getViz(mod8)) +
         ggtitle("s(Distance) + s(Windspeed) + s(Precipitation)"),
      pages = 1)

0.2.9 mod9 - s(Distance) + s(Precipitation)

mod9 <-
   gam(
      mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5),
      data = dat
   )

summary(mod9)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   1.0802     0.0459    23.5 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F              p-value    
## s(distance) 3.93    4.0 83.8 < 0.0000000000000002 ***
## s(sum_rain) 2.03    2.2 10.8             0.000023 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.515   Deviance explained = 52.4%
## GCV = 0.71995  Scale est. = 0.70494   n = 334
print(p_gam(x = getViz(mod9)) +
         ggtitle("s(Distance) + s(Precipitation)"),
      pages = 1)

0.2.10 mod10 - s(Distance) +s(Precipitation) + Windspeed

mod10 <-
   gam(
      mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5) + mws,
      data = dat
   )

summary(mod10)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5) + mws
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -2.863      1.464   -1.96   0.0513 . 
## mws            1.109      0.412    2.70   0.0074 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df    F              p-value    
## s(distance) 3.94   4.00 93.8 < 0.0000000000000002 ***
## s(sum_rain) 3.82   3.97 11.2         0.0000000097 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.566   Deviance explained = 57.8%
## GCV = 0.64947  Scale est. = 0.6305    n = 334
print(p_gam(x = getViz(mod10)) +
         ggtitle("s(Distance) + s(Precipitation) + Windspeed"),
      pages = 1)

0.2.11 mod11.0 - s(Distance) + s(Windspeed) + s(Precipitation), family = tw()

This is the same as mod8 but using family = tw(), see ?family.mgcv for more on the families. The Tweedie distribution is used where the distribution has a positive mass at zero, but is continuous unlike the Poisson distribution that requires count data. The data visualisation shows clearly that the mean pot count data have this shape.

mod11.0 <-
   gam(
      mean_pot_count ~ s(distance, k = 5) + 
         s(mws, k = 5) + 
         s(sum_rain, k = 5),
      data = dat,
      family = tw()
   )

summary(mod11.0)
## 
## Family: Tweedie(p=1.044) 
## Link function: log 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5) + s(sum_rain, 
##     k = 5)
## 
## Parametric coefficients:
##             Estimate Std. Error t value    Pr(>|t|)    
## (Intercept)   -0.228      0.041   -5.57 0.000000053 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df      F              p-value    
## s(distance) 3.50   3.85 123.76 < 0.0000000000000002 ***
## s(mws)      2.94   3.05  13.50             0.000014 ***
## s(sum_rain) 1.90   1.93   5.57                0.013 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.673   Deviance explained = 61.2%
## -REML = 309.77  Scale est. = 0.36396   n = 334
print(p_gam(x = getViz(mod11.0)) +
   ggtitle("s(Distance) + s(Windspeed) + s(Precipitation), family = tw()"),
   pages = 1)

0.2.12 mod11.1 - s(Distance, bs = “ts”) + s(Windspeed, bs = “ts”) + s(Precipitation, bs = “ts”), family = tw()

mod11.1 <-
   gam(
      mean_pot_count ~ s(distance, k = 5, bs = "ts") + 
         s(mws, k = 5, bs = "ts") + 
         s(sum_rain, k = 5, bs = "ts"),
      data = dat,
      family = tw()
   )

summary(mod11.1)
## 
## Family: Tweedie(p=1.044) 
## Link function: log 
## 
## Formula:
## mean_pot_count ~ s(distance, k = 5, bs = "ts") + s(mws, k = 5, 
##     bs = "ts") + s(sum_rain, k = 5, bs = "ts")
## 
## Parametric coefficients:
##             Estimate Std. Error t value   Pr(>|t|)    
## (Intercept)  -0.2259     0.0409   -5.52 0.00000007 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##              edf Ref.df     F              p-value    
## s(distance) 3.25      4 117.9 < 0.0000000000000002 ***
## s(mws)      2.88      4  13.4      0.0000000000019 ***
## s(sum_rain) 1.84      4   3.0               0.0005 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.671   Deviance explained =   61%
## -REML = 322.33  Scale est. = 0.36412   n = 334
print(
   p_gam(x = getViz(mod11.1)) +
      ggtitle(
         "s(Distance, bs = 'ts') + s(Windspeed, bs = 'ts')\n+ s(Precipitation, bs = 'ts'), family = tw()"
      ),
   pages = 1
)

This model, same structure as mod11.0, uses thin-plate splines to shrink the coefficients of the smooth to zero when possible.

0.3 Compare the Models

0.3.1 AIC, BIC

models <- list(mod1 = mod1,
               mod2 = mod2,
               mod3 = mod3,
               mod4 = mod4,
               mod5 = mod5,
               mod6 = mod6,
               mod7 = mod7,
               mod8 = mod8,
               mod9 = mod9,
               mod10 = mod10,
               mod11.0 = mod11.0,
               mod11.1 = mod11.1
               )
map_df(models, glance, .id = "model") %>%
   arrange(AIC)
## # A tibble: 12 x 7
##    model      df logLik   AIC   BIC deviance df.residual
##    <chr>   <dbl>  <dbl> <dbl> <dbl>    <dbl>       <dbl>
##  1 mod11.0  9.34  -320.  663.  708.     141.        325.
##  2 mod11.1  8.96  -322.  667.  712.     141.        325.
##  3 mod10    9.76  -392.  805.  846.     204.        324.
##  4 mod8     9.85  -392.  805.  847.     204.        324.
##  5 mod6     9.87  -392.  806.  847.     204.        324.
##  6 mod7     8.85  -394.  808.  845.     207.        325.
##  7 mod9     6.96  -412.  840.  870.     231.        327.
##  8 mod4     6.93  -414.  844.  874.     233.        327.
##  9 mod3     5.93  -416.  846.  873.     236.        328.
## 10 mod2     5.93  -423.  860.  886.     246.        328.
## 11 mod5     5.93  -423.  860.  886.     246.        328.
## 12 mod1     4.93  -424.  860.  883.     248.        329.

0.3.2 R2

enframe(c(
   mod1 = summary(mod1)$r.sq,
   mod2 = summary(mod2)$r.sq,
   mod3 = summary(mod3)$r.sq,
   mod4 = summary(mod4)$r.sq,
   mod5 = summary(mod5)$r.sq,
   mod6 = summary(mod6)$r.sq,
   mod7 = summary(mod7)$r.sq,
   mod8 = summary(mod8)$r.sq,
   mod9 = summary(mod9)$r.sq,
   mod10 = summary(mod10)$r.sq,
   mod11.0 = summary(mod11.0)$r.sq,
   mod11.1 = summary(mod11.1)$r.sq
)) %>%
   arrange(desc(value))
## # A tibble: 12 x 2
##    name    value
##    <chr>   <dbl>
##  1 mod11.0 0.673
##  2 mod11.1 0.671
##  3 mod10   0.566
##  4 mod6    0.566
##  5 mod8    0.566
##  6 mod7    0.562
##  7 mod9    0.515
##  8 mod4    0.509
##  9 mod3    0.504
## 10 mod2    0.484
## 11 mod5    0.484
## 12 mod1    0.482

0.3.3 ANOVA

anova(mod1,
      mod2,
      mod3,
      mod4,
      mod5,
      mod6,
      mod7,
      mod8,
      mod9,
      mod10,
      mod11.0,
      mod11.1,
      test = "F")
## Analysis of Deviance Table
## 
## Model  1: mean_pot_count ~ s(distance, k = 5)
## Model  2: mean_pot_count ~ sum_rain + s(distance, k = 5)
## Model  3: mean_pot_count ~ mws + s(distance, k = 5)
## Model  4: mean_pot_count ~ sum_rain + mws + s(distance, k = 5)
## Model  5: mean_pot_count ~ sum_rain + s(distance + mws, k = 5)
## Model  6: mean_pot_count ~ sum_rain + s(distance, k = 5) + s(mws, k = 5)
## Model  7: mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5)
## Model  8: mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5) + s(sum_rain, 
##     k = 5)
## Model  9: mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5)
## Model 10: mean_pot_count ~ s(distance, k = 5) + s(sum_rain, k = 5) + mws
## Model 11: mean_pot_count ~ s(distance, k = 5) + s(mws, k = 5) + s(sum_rain, 
##     k = 5)
## Model 12: mean_pot_count ~ s(distance, k = 5, bs = "ts") + s(mws, k = 5, 
##     bs = "ts") + s(sum_rain, k = 5, bs = "ts")
##    Resid. Df Resid. Dev        Df Deviance        F               Pr(>F)    
## 1        329        248                                                     
## 2        328        246  1.000045        2     4.84               0.0285 *  
## 3        328        236  0.000324       10 83384.02    0.000000000027800 ***
## 4        327        233  1.000093        3     8.75               0.0033 ** 
## 5        328        246 -1.000418      -13    35.79    0.000000005797152 ***
## 6        324        204  3.997805       42    28.85 < 0.0000000000000002 ***
## 7        325        207 -1.001943       -3     7.17               0.0077 ** 
## 8        324        204  0.998301        3     7.16               0.0079 ** 
## 9        327        231 -2.795499      -26    25.71    0.000000000000035 ***
## 10       324        204  2.769762       26    25.90    0.000000000000036 ***
## 11       324        640  0.372293     -435                                  
## 12       324        644 -0.015787       -5   811.78    0.000002140998087 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.3.4 Check Best Model Fit

0.3.4.1 Mod11.0 - s(Distance) + s(Windspeed) + s(Precipitation), family = tw()

mod11.0_vis <- getViz(mod11.0)
check(mod11.0_vis,
      a.qq = list(method = "tnorm", 
                  a.cipoly = list(fill = "light blue")), 
      a.respoi = list(size = 0.5), 
      a.hist = list(bins = 10))
## 
## Method: REML   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-0.0006677,-0.000000007453]
## (score 309.8 & scale 0.364).
## Hessian positive definite, eigenvalue range [0.394,2978].
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value   
## s(distance) 4.00 3.50    0.86   0.010 **
## s(mws)      4.00 2.94    0.89   0.045 * 
## s(sum_rain) 4.00 1.90    0.90   0.045 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.3.4.2 Mod11.0 - s(Distance, family = “ts”) + s(Windspeed, family = “ts”) + s(Precipitation, family = “ts”), family = tw()

mod11.1_vis <- getViz(mod11.1)
check(mod11.1_vis,
      a.qq = list(method = "tnorm", 
                  a.cipoly = list(fill = "light blue")), 
      a.respoi = list(size = 0.5), 
      a.hist = list(bins = 10))
## 
## Method: REML   Optimizer: outer newton
## full convergence after 8 iterations.
## Gradient range [-0.000000006816,0.0000000575]
## (score 322.3 & scale 0.3641).
## Hessian positive definite, eigenvalue range [0.6649,2980].
## Model rank =  13 / 13 
## 
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
## 
##               k'  edf k-index p-value   
## s(distance) 4.00 3.25    0.86   0.010 **
## s(mws)      4.00 2.88    0.88   0.020 * 
## s(sum_rain) 4.00 1.84    0.89   0.025 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

0.4 Thoughts

This model, mod11.0, mean_pot_count ~ s(Distance) + s(Windspeed) + s(Precipitation) - family = tw(), is the best performing model. It cannot be used for predictions, but suitably describes the dispersal data we have on hand with the parameters used. More data would be desireable to increase the value of k as evidenced in the GAM checks.